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Abstract. Glauber dynamics of a bond-diluted Ising model on a Bethe lattice (a 
random graph with fixed connectivity) is investigated by an approximate theory which 
provides exact results for equilibrium properties. The time-dependent solutions of 
the dynamical system derived by this method are in good agreement with the results 
obtained by Monte Carlo simulations in almost all situations. Furthermore, the derived 
dynamical system exhibits a remarkable phenomenon that the magnetization shows 
multi-step relaxations at intermediate time scales in a low-temperature part of the 
Griffiths phase without bond percolation clusters. 
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1. Introduction 

Clarifying the role of impurities in many-body systems is a central issue in statistical 
physics because statistical quantities of a system with a kind of impurities are often 
different from those of a system without the impurities qualitatively. Griffiths-McCoy 
singularities, corresponding to the essential singularity appearing in the density of Lee- 
Yang zeros of a partition function, are such typical phenomena [1, 2, 3]. 

Thus far, the equilibrium properties of Griffiths-McCoy singularities have been 
extensively investigated [3, 4]. In order to simplify the problem on finite-dimensional 
systems, models on Bethe lattices have been studied with special theoretical techniques, 
which accelerate the understanding of the equilibrium properties of the singularities 
[5, 6, 7, 8]. However, in contrast to the considerable amount of knowledge available 
on equilibrium properties, the knowledge available on dynamical properties associated 
with Griffiths-McCoy singularities is scarce except for a few previous studies which have 
clarified anomalous dynamical behaviours in the long-time limit [9, 10, 11]. Thus, the 
knowledge of dynamical behaviours at intermediate time scales has not been established 
although such a challenge for a model on a Bethe lattice exists [12]. However, even 
if a Bethe lattice is adopted as a simple situation, the understanding of dynamical 
behaviours at intermediate time scales is still limited. Probably, one of the reasons 
for this limited understanding is the lack of theoretical treatments for such dynamical 
problems. 

In this study, in order to solve such a problem, we attempt to develop an 
approximate theory for analysing the dynamics of a bond-diluted Ising model on a 
regular random graph with fixed connectivity (a Bethe lattice). We find that the final 
states of the derived dynamical system are exact, and the time-dependent solutions of the 
dynamical system are in good agreement with the results obtained by Monte Carlo (MC) 
simulations at almost all situations, although there are slight discrepancies between 
the two results below the Griffiths-paramagnet transition temperature with bond 
percolation clusters. From the results of this method, we predict that magnetization 
shows multi-step relaxations at intermediate time scales in a low-temperature part of 
the Griffiths phase without bond percolation clusters. 

This paper is organized as follows. In section 2, we define a bond-diluted Ising model 
precisely. In section 3, we derive a dynamical system from the model approximately. In 
section 4, we investigate the concrete behaviours of the derived dynamical system and 
compare them with MC simulations. In the concluding remarks, we discuss the validity 
of the approximate method. In the Appendix, we present the equilibrium properties of 
the bond-diluted Ising model on the Bethe lattice. 

2. Model 

Let us consider a regular random graph, which consists of N G N sites, and each 
site connects to c e N sites chosen randomly where N is the set of natural numbers. 
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Then, G(c,N) is denned as a set of such regular random graphs. For the spin variable 
<jj G {—1, 1} defined on each site % G {1, • • • , N} in a random graph Q G G(c, N), we 
consider a bond-diluted Ising model, whose Hamiltonian is given as 

^V) = -^EE J wi> (!) 

where 5j is defined as a set of sites connected to site i, and we collectively express 
a = [?i)i=\ an d J = (Jij)fj =1 . The configuration of bonds (= J,-j) is given by the 
probability 

D(Jn) = P 5(J tJ - 1) + (1 - p)S(Jij), (2) 

for j G Sj, otherwise Jij = 0. The average of a quantity X over diluted bonds 
is expressed as X. Next, we define the dynamics of the bond-diluted Ising model 
as follows. Let T J {Fi<r — \ cr) be the transition rate from i^cr to cr which satisfies 
the detailed balance condition, where Fi is the spin flip operator at site % such that 
Fi<r = (<7i, • • • , — <7j, • • • , (Tat). The master equation for the probability P J (cr, t), that 
the spin configuration is cr at time t with a realization of J, is as follows: 

N 

d t P J ((T,t) = Y,[T J (F l( T Cr)P J (F,CT,t) - T J (cr -> F t Cr)P J (cr, *)]. (3) 
i=l 

In this paper, we consider the case 

T J (F t a cr) = - tanh((if J (F l cr) - // J (cr))/2T)]. (4) 

Here we briefly review the equilibrium properties of Griffith-McCoy singularities. 
Let Z{0) be a partition function of a pure imaginary external field 6. The term of 
Griffiths-McCoy singularities is defined as the fact that in the thermodynamic limit, 
the density of zeros of Z(6) in the imaginary axis has the essential singularity such as 
exp(—A(p,T)/\9\) in the limit \9\ — > 0, where A(p,T) is a constant depending on the 
values of parameters (p, T). It should be noted that with the values of the parameters at 
which such singularities occur, the magnetization in equilibrium with 9 = is still zero. 
The Griffiths phase where such singularities are observed is shown in figure 1. The phase 
boundary in terms of the ferromagnetic phase is obtained from the discussion presented 
in Appendix A. It should be noted that the phase boundary between the Griffiths phase 
and the paramagnetic phase is not obtained by the discussion presented in Appendix A. 
More details about Griffiths- McCoy singularities are discussed in [8]. Naively speaking, 
Griffiths-McCoy singularities originate from rare large ordered regions. The main focus 
of this study is to understand the anomalous dynamical behaviours associated with the 
singularities. In the next section, we construct an approximate theory to understand 
such behaviours. 

3. Derivation of an effective dynamical system 

In principle, master equation (3) provides complete information about the system. 
However, it is difficult to extract useful information from this equation because the 
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Figure 1. Phase diagram of the system with c = 3 where F, P and G indicate the 
ferromagnetic phase, paramagnetic phase and Griffiths phase, respectively. 



number of states of the system is 2 , which is quite large when N is large. Furthermore, 
the analysis of the system becomes more complicated because of the presence of diluted 
bonds as a quenched disorder. Here, let us remind that a useful theoretical method 
has been constructed to understand equilibrium properties of some systems on Bethe 
lattices without directly considering such large number of variables [13]. On the basis of 
the effectiveness of this method, first, we attempt to describe the system with a set of a 
finite number N c g G O(N ) variables in order to obtain more useful information on the 
dynamical behaviours of the system, However, in general, it is not straightforward to 
extend the exact analysis from equilibrium properties to dynamical properties. In fact, 
although some works have been succeeded in deriving the exact evolution equations for a 
reduced number of variables [14, 15, 16, 17, 18, 19], each method has each special aspect, 
by which it is difficult to obtain useful information about the dynamical properties of 
the system in the case of the present problem. As another way to obtain them, let us 
focus on a remarkable previous study [20]. The study has reported that an approximate 
description with a finite number of variables, which is exact for equilibrium states, can 
be obtained in an Ising ferromagnet. Here, on the basis of the results of this study, we 
attempt to develop an approximate theory for dealing with the quenched disorders of 
the system on Bethe lattices. 

First, as a preliminary step, we call site j G Bi with the bond = 1 'bonded site 
of site i\ and we call site j G Bi with the bond = 'unbonded site of site V . The 
number G {0, ■ ■ • , c} of bonded sites of site i is determined by 



Let Ui G {0, ■ ■ • , U} denote the number of downward spins at the bonded site of site i. 
Then, site i is characterized by a set of variables (<7j, Ui, k). Using these variables, we can 
express T J (cr — > Fjcr) as r ai (li,Ui), where r a (l,u) = (1— tanh(cr(Z — 2u)/T))/2. It should 
be noted that <7j and depend on time t, whereas k does not depend on time t. Next, 
with a time t and a realization of diluted bonds «J, let P/(o~, u; t) be the probability that 
(crj,Ui) takes (a, u), and let PJ((a,u), (cr',u');t) be the joint probability that (o"j,Uj) 




(5) 



Glauber dynamics of a bond- diluted Ising model on a Bethe lattice 5 

and ((Tj,Uj) take (a,u) and (a', u'), respectively. Using these probabilities, we also 
define the conditional probability P^(a,u\a' \u';t) = Pf j ((a,u),(a',u');t)/Pj r (a',u';t). 
In the following analysis, we fix a graph Q e G(c,N) for sufficiently large N without 
considering the ensemble for G(c,N). That is, the following analysis can be applicable 
to almost all graph Q e G(c,N) in the thermodynamic limit. 

With these notations, we start with the following exact expression: 

d t Pf(a, u- 1) = -r a (u, h)Pf(a, u; t) + r_ a («, U)Pf{-a, u; t) 

h 

+ r a ,(u j ,l j )Pf(a',u j ;t)6(J ij - 1) 

j£Bi a' Uj=0 

[P*{a, u + a'\a', Uj] t) - Pf^a, u\a', Uj ; t)}, (6) 

where we used the property that the order of the length of loops in the random graph 
Q is O(logiV) longer than three, and we modified the definition P J (a, u\a', c + l;t) = 
P J (a,u\a', — l;t) = for convenience. We carry out an average over diluted bonds by 
multiplying both the sides of equation (6) with 5(^2j eB . Jij — I). Here, we define 



Pi(<r, I, u; t) = 5(J2 Jij - O^V, «; *), (7) 

jeBi 

Pij(o~, I, u\a', I', u'\ t) = 

5(^-1)5(^2 Jik-l)8( Yl Jjk'-l'^o-Av'iW-it), (8) 



Pj(a, I, u; t) = 6(Jii - 1)5(J2 Jij ~ l)P?{o, «l *)■ (9) 

jeBi 

Then, we have 

d t Pi((T, I, u; t) = -r a (l, u)Pi(a, I, u; t) + r_ CT (Z, u)Pi(-a, I, u; t) 

jeBi a' u'<V 

[Pij(a, l,u + a'\a', I', u'\ t) - P^a, I, u\a', I', u; t)}, (10) 

where ^2 ul<l , = Y^'=o Y^ u '=cr identities for the derivation of (10) from equation (6) 
are 



c 



^ = E ,5 (E j «- z ')a, (ii) 

l'=o jeBi 

-5(a)Y(a) 



5(a)X(a)Y(a)/Z(a) = 5(a)X(a)=^=^. (12) 

o(a)Z{a) 

Hereafter, we assume that the values of P](o' ,1' ,u';t) and Pij(a,l,u'\a',l',u';t) do 
not depend on the two chosen bonded sites i,j in the thermodynamic limit. This 
assumption may be plausible because at least, in this case, inhomogeneous properties of 
the system originated from the effects of loops of the random graph Q may be negligible 
in the thermodynamic limit. This assumption corresponds to that Pij(cr, /, u'\a', I', u'\ t) 
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is the same as the conditional probability Pb{&, I, u\a' ', I', u'; t) that if a site characterized 
by (<t, /, u) is randomly chosen, another site characterized by (a, I, u) is randomly chosen 
in the bonded sites of the first chosen site, and also that Pj(o~', I', u'; t) is equal to 
{V /c)Pj(o J ,1' ,u';t). Here, we define P^(l,u\a',l',u';t) as follows: 

TO,|^,.';«) S ^y f ^;" , (13) 

Pb{o-\u ', L',U , t) 

where Pb(c|c', l', u '\ t) is defined by the similar way to that of Pb(c, I, u'\a', I', u'\ t). 
Then, we can obtain the following relation: 

P B («rK, f , t) = Mf^, (14) 

°*<* [:=-\, < i6 > 

In the above procedures with a trivial relation (52j eB .) = c, dynamical system (10) is 
simplified into 

d t Po,i,u{t) = -r a (l,u)p a ^ u (t) + r_ <T (Z,u)p_ ff>J)tt (*) 

+ ^^v(!» fl/y (f)a(f,«') 

a' u'<V 

[P°{l,u + oy,l\u'-t)-P°{l,u\o-',l',u'-t)\, (16) 

where p a , u ,i(t) = ^2iLiPi(o~,u,l;t)/N. It may be plausible that p a , u ,i{t) is identical to 
Sili ~~ °"i(^))^( M — Ui(t))S(l — k)/N in the thermodynamic limit N — > oo. At the 
present level, equation (16) is not closed in terms of p = (p a ,i, u )- 

Next, in order to obtain a closed equation based on (16), we carry out the following 
approximation in equation (16): 

P£(l, u\a', V, u'\ t) = P°(l, u\a'] t). (17) 

By considering the number of connected sites on which the spin variable is a', we can 
express Pb,(t(1, u\a'; t) using p as follows: 

m uW t) = h l ~ u )p*jM/ <r p**At) f = +i) (18) 
Bl 1 ; |«toW/Ev<,'4,iv(«) K = -i)- 1 ; 

From (16) with approximation (17), we obtain a closed evolution equation in terms of 
p as follows: 

dtP = G(p), (19) 

where G = (G a j tU ), and G a j >u is a function expressed by the right-hand side of equation 
(16) with approximation (17) and relation (18). At the present level, the number of 
system variables is reduced to N eS = (c + l)(c + 2)/2 from 2^. 

In order to consider the validity of approximation (17), first, we focus on the 
stationary solution p cq satisfying G(p eq ) = and define the time t eq = mint st such 
that p(t st ) = p cq . By rewriting 

P°(l uW V u'- 1) - PB({o-',l',u'),{o-,l,u);t) 
B (/ ' Ul ' 1 ' U ' } " £„„<,* Pb((o->, l\ u>), (a, I", «"); t) ' (2 ° } 
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where 

P B ((a,l,u),(a',l',u');t)= 

5(J tJ - 1)S(J2 Jik ~ l)S( J i« ~ (a',u');t), (21) 

with the assumption that the value of Pb((o~, l,u), (a',l f ,u');t) does not depend on the 
chosen bonded sites we find that in the equilibrium case t > t cq , Pg(l, u\a', I', u'\ t) = 
P£(l,u\a';t) holds because of relations (18), (B.6) and (B.7). (See Appendix B for the 
details of the derivation of relations (B.6) and (B.7)) That is, in the equilibrium case 
t > t cq , approximation (17) becomes an identity. 

It should be noted that if cr at t = is given, p(0) is automatically determined. 
That is, in the case m(0) = 1, p(0) is given as 

P+i,i,o(0) = Q jp'(l-p) c -', (22) 

P-Uo(O) =0, (23) 

P*,i,uM$) = 0- (24) 

We mention that the magnetization m = o~i/N and the energy density e = H(&)/N 
can be expressed in terms of p as follows: 

a u<l 

e =-^$>(Z-2 M ) P(tM . (26) 

(J U<1 

A conjecture arising from the used approximation 

We consider the properties of used approximation (17) for the relaxation processes 
with the initial condition m(0) = 1. Here, let m and r be the magnetization and the 
relaxation time of the model in the thermodynamic limit. Further, let and Td be the 
magnetization and the relaxation time described by the solutions of dynamical system 
(19). Within the used approximation, we ignore the correlations of two sites between 
which length is far than two sites, although such correlations have finite relaxation times. 
This leads to the conjecture that there exists a lower bound for the relaxation time r 
such that rd < r, which might lead to a lower bound m&it) < m(t) for arbitrary t with 
m(0) = 1. Partially, this conjecture is supported by the systematic work considering 
the effects of far sites in an Ising model with a ferromagnetic coupling [20]. In future, 
it will be important to provide mathematical proof of this conjecture. 

4. Time-dependent solutions of dynamical system (19) 



We analyse the time-dependent solutions of the derived dynamical system in detail. We 
consider the relaxation behaviours of the system from the initial condition m(0) = 1 
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Figure 2. Time-dependent magnetization described by dynamical system (19) and 
MC simulation from the initial condition m(0) = 1. (left) p = 1 (right) p = 0.7. In the 
figure, MC denotes MC simulations, D denotes dynamical system (19). 



with fixed values of parameters (p,T). Concretely, we numerically solve dynamical 
system (19) using the forth-order Runge-Kutta method with a time step St = 10~ 2 . In 
this paper, we mainly consider the system with c = 3. In fact, the system with c > 3 
exhibits the same behaviours as the system with c = 3 qualitatively, and the system 
with c = 2 is equivalent to a one-dimensional system in which any phase transitions do 
not occur. 

This section is organized as follows. In section 4.1, we focus on the Ising ferromagnet 
limit p — 1, and we confirm that the results are consistent with those of a previous 
work. In section 4.2, we focus on the ferromagnet phase and the paramagnet phase with 
p > p c . In section 4.3, we focus on the Griffiths phase with p > p c , and we find a slight 
discrepancy between the dynamical system and the MC simulations. In section 4.4, we 
focus on the Griffiths phase with p < p c , and we find noteworthy behaviours associated 
with magnetization. It should be noted that MC simulations are carried out by applying 
the following dynamical rule. First, a site i is randomly chosen. Then, the spin variable 
0-j is flipped at the probability T J (<x — y Fi<r). This step is repeated, and time t = 1 is 
defined for such repeated iV steps. This dynamical rule is equivalent to the dynamics 
of the model that we consider in the thermodynamic limit. In all MC simulations, the 
system size N is 10 6 , the number of samples for the statistical average is 10 where each 
sample has a different realization of J and G(c, N). 

4-1. Ising ferromagnet limit 

Dynamical system (19) with p = 1 exactly corresponds to the dynamical system 
derived by the 'Independent-neighbor approximation' proposed in the previous work 
[20]. According to the work, this dynamical system can capture critical exponents of 
the dynamical critical phenomena at Tq(1). For example, the dynamical system can 
provide the relaxation time r of the system exhibiting r ~ e± with £ = 1 where 
e± = ±(T — Tq(p)). In the left-hand side of figure 2, we show two examples of the time- 
dependent magnetization described by (19) together with the results of MC simulations. 
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Figure 3. Time-dependent magnetization described by dynamical system (19) in the 
Griffiths phase with T — 1.2, 1.3 from the initial condition m(0) = 1. c = 3, p — 0.7. 
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Figure 4. Time-dependent magnetization described by dynamical system (19) and 
MC simulations in the Griffiths phase with T = 1.2,1.3 from the initial condition 
m(0) = 1. c = 3, p = 0.7. 



It is seen that the coincidences between (19) and MC simulations are quite good. 
4-2. The ferromagnetic phase, and the paramagnetic phase 

Here, we focus on the behaviours with the parameter p = 0.7. In the right-hand side of 
figure 2, we show two examples of the time-dependent magnetization described by (19) 
together with the results of MC simulations. It is seen that the coincidences between 
(19) and MC simulations are quite good, although there is a slight discrepancy between 
two results in the ferromagnetic phase T = 0.9. In the concluding remarks, we discuss 
the discrepancy. 

4-3. The Griffiths phase with bond percolation clusters 

We focus on the behaviours with the parameter p = 0.7. As shown in figure 3, the 
relaxational form of the magnetization seems to be nonexponential as a function of time 
at an intermediate time scale; however it exhibits an exponential in the long-time limit. 
In fact, the signs for nonexponential decay are found by MC simulations as shown in 
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Figure 5. Time-dependent magnetization described by dynamical system (19) from 
the initial condition m(0) = 1 with c = 3, p = 1/5. 
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Figure 6. Each quantity described by dynamical system (19) from the initial condition 
to = 1 at T = 0.4 with c = 3, p = 1/5. 



figure 4. This means that the used approximation by the analysis in this paper fails to 
capture such anomalous behaviours in the long-time limit, although the approximation 
works well until an intermediate time scale. In the concluding remarks, we discuss the 
validity of the used approximation. It should be noted that the existence of such a 
slow relaxation has been rigorously proved in the models on <i-dimensional lattices [10] 
although as far as we know, there is no rigorous proof of such an argument for the 
models on Bethe lattices. 



4-4- The Griffiths phase without bond percolation clusters 

Next, we focus on the case p = 0.2. Near the transition point Tq(1), the relaxation 
form of the magnetization is exponential in a wide region, as seen in the left-hand side 
of figure 5. When the temperature is decreased to around T = 0.6, contrary to such 
behaviours, the magnetization surprisingly shows the multi-step relaxation as seen in 
the right-hand side of figure 5. Rigorously, we use the term '/-step' when there exists 
the time tj for j e {1, ■ ■ • ,/} such that <9 2 (logm)/<9(logt) 2 |t = f. = 0. In this sense, the 
three-step relaxation occurs in the case c = 3. 

We expect that this observation can be related to the existence of the following 
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Figure 7. Time-dependent magnetization described by dynamical system (19) from 
the initial condition m(0) = 1 with c = 5,p = 1/10. 

terms in (19); 

9tP+i,o-ifi(t) = - r +i( c - *: 0)p + i, c _i, (t) + A(p(t)), (27) 
where i G {0, ■ ■ • , c}, and A(p(t)) is a function of p(t). That is, each variable has different 
relaxation time Tj ~ r + i(c— i, 0) _1 , which is almost exp((c— i)/T) with T< 1, assuming 
that the dependence of on p + i )C _ i)0 (£) can be ignored. In fact, as shown in figure 

6, we confirm that each relaxation time r, is very close to each inflection point defined 
above on the relaxation form of the magnetization. Our conjecture is that c + 1 number 
of the relaxation times lead to c number of the inflection points. For the case c = 4,5, 6, 
we have confirmed that the c-step relaxation occurs in dynamical system (19). In the 
figure 7, we show the five-step relaxation in the magnetization for the case c = 5. 

However, with only this consideration, we cannot understand the mechanism of the 
multi-step relaxations in the magnetization because we have not yet understood what 
determines the threshold temperature 2th where a multi-step relaxation begins to occur. 
We probably believe that a key fact for the appearance of the multi-step relaxations is 
that the relative rate of each quantity, for example, p+i,i,i(ieq)> P+i,2,i(£eq), P+i,3,o(£eq) as 
a function of temperature T is qualitatively changed around T th , as shown in figure 
8. Concretely, when the temperature is decreased, p+i,3,o(£eq) is increased in the 
manner p +ljljl (i(; cq ) > p +1>2 ,i(t cq ) > p+i,3,o(*eq) in high-temperature regions, p+i,i,i(t e q) > 

P+l,3,o(*eq) > P+l,2,l(*eq) around T th ~ 0.6, and p+l, 3 ,o(*eq) > P+l,l,l(*eq) > P+l,2,l(*eq) in 

low temperature regions. This leads to the conjecture that this qualitative change in the 
relative rates drives each relaxation time Tj become apparent in the magnetization as a 
multi-step relaxation. For the case c = 4,5, 6, we have confirmed that similar changes 
in the inequality to the case c = 3 occur around T th (p, c) depending on the parameters 
p and c. In spite of such observations, we still cannot understand how the change in the 
inequality is related to the appearance of the multi-step relaxations; in future, we intend 
to conduct a study to clarify it. It might also be an interesting to clarify the mechanism 
of the multi-step relaxations by performing a theoretical analysis of dynamical system 
(19) directly. 

In figure 9, we show two examples of the time-dependent magnetization described 



Glauber dynamics of a bond- diluted Ising model on a Bethe lattice 



12 






Figure 8. p+i,i,i(t e q), P+i,2,i(ieq) and jO+i,3,o(*eq) depending on temperature T with 
c = 3, p = 1/5. 
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Figure 9. Time-dependent magnetization described by dynamical system (19) and 
MC simulation in Griffiths phase from the initial condition m(0) = 1. c = 3, p — 1/5. 



by (19) together with the results of MC simulations. It is seen that the coincidences 
between (19) and MC simulations are also quite good even if a multi-step relaxation 
occurs. This result strongly suggests that the 'true' magnetization also shows a multi- 
step relaxation even in the case c = 5 such as the magnetization described by dynamical 
system (19) in figure 7. In this parameter, it is difficult to judge by MC simulations 
whether such behaviors occur or not. Finally, we mention that even if such multi-step 
relaxations occur, the relaxation form of the magnetization in the long-time limit is 
expected to be exponential on the basis of the form (27). It should be noted that the 
existence of such an exponential form has been rigorously proved in the models on d- 
dimensional lattices [11] although as far as we know, there is no rigorous proof of such 
an argument for the models on Bethe lattices. 

5. Concluding remarks 

In this section, we mention the relevance of our study to the other studies. So far, as an 
extension of the previous work [20], a method called dynamical replica theory (DRT) 
has been already developed to analyse the system on Bethe lattices with quenched 
disorders [12, 21]. For the bond-diluted Ising model, the dynamical replica theory has 
already clarified the existence of a special temperature, where the relaxation of the 
system is qualitatively changed, and the qualitative change in the relative rate of some 
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spin-field distributions [12]. However, the concrete form of the relaxation below the 
special temperature has not been clear yet in the study [12]. The analysis in this paper 
has proceeded the understanding of this phenomenon. That is, it has clarified that 
the magnetization below the special temperature shows multi-step relaxations. It is an 
important future study to clarify the relation between the method in this paper and the 
method of DRT because the two method may be complementary at each other. 

Before concluding this paper, we discuss the validity of the used approximation. 
We mention that the used approximation also fails to capture the anomalous 
relaxation dynamics of a kinetically constrained Ising model on a Bethe lattice, where 
heterogeneous spin configurations ordered by a kinetic constraint exist [22]. Therefore, 
our conjecture is that the failure in capturing the anomalous relaxation behaviours 
below the Griffiths-paramagnetic transition temperature with bond percolation clusters 
is originated by some heterogeneous spin configurations in the relaxation dynamics, 
which might be related to rare large ordered regions. Furthermore, we have evidences 
for the conjectures that we can obtain better dynamical systems gradually for some 
systems by enlarging the environmental regions of a site, which are considered by 
effective variables [20, 22]. Therefore, we expect that in principle, such a procedure can 
gradually extend the time span where nonexponential behaviours are observed as found 
in section 4.3, although such analysis is very complicated. Finally, we mention that 
it is an open problem how the used assumption related to inhomogeneous properties 
originated from the graph structure and the used approximation work at the system 
where the replica symmetry breaking occurs. 
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Appendix A. The critical temperature Tq{p) 

Let us consider a Cayley tree, which has the same structures as those of a random 
graph Q G G(c,N), ignoring the loops of which length is O(logiV). Here, we denote 
the generation of the Cayley tree by g G {0, 1, • • • , g max }, where g — is assigned 
to the root. For a given site in (g — l)-th generation, a bond connecting the site 
with a site in g-th generation is labeled by an integer j g G {1, • • • , c — 1}. Thus, an 
arbitrary site in g-th generation of the Cayley tree is indicated by a set of integers 
j g = U 9 ,---,ji,0). Let C(j g ) be a set of sites defined as {(i 9majc , ■ ■ ■ , i g+1 ,j g )\l < 
ik < c - l,g + 1 < k < g max } where (i,k,j g ) = (i, k, j g , ■ ■ ■ , j 1 , 0). The sites 
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in C(j g ) can be regarded as sites on a Cayley tree T 3 - 9 1 with the root being site 
j g . Then, for a realization of diluted-bonds J, let us consider a partition function 
Z J jg {°) = E ie cy 9 )E^K ~ ^expi-PH'iWiecuj)) where B = l/T. Here, we 
define exp(/32/if ) = (— 1)/Z? (1). Then, we can obtain the recursive equation in 

■> 9 ■> 9 3 _. 

terms of the distribution Qj g (h) = S(h — h J - ) as follows: 

Q jg (h)=H([Q jg+l (h)}), (A.l) 

H([Q(h)}) = / dhiQihiMh-YtUVwj.'hi), (A.2) 
J i=i 

where u(J,h) = Q~ x tanh _1 [tanh(/3 J) tanh(/3/i)]. 

It should be noted that by solving recursive equation (A.l) with almost all values of 
initial conditions Qj g (h) with g = g max , Qj (h) for arbitrary j g with g <C g max becomes 
a solution Q(h) satisfying the following equation 

Q(h) = H([Q(h)}). (A.3) 

According to the previous study [8], for a sufficiently high temperature, the solution Q(h) 
is 5(h). When the temperature T is decreased from a sufficiently high temperature 
for a fixed p, we can obtain the critical temperature Tq(p) where the point hj g (= 
f dhQj (h)h) = loses linear stability in the recursive equation between hj g and hj , 
which is derived from recursive equation (A.l). From this, we obtain the relation 
1/Tq(p) = tanh _1 (p c /p), where p c = l/(c— 1) is a threshold value in a bond-percolation 
problem. That is, there are bond percolation clusters in the system with p > p c , no 
bond percolation clusters in the system with p < p c . Concretely, Tq(1) = 1.820478 for 
c = 3, and for c < 3, any phase transitions do not occur in the model. 



Appendix B. Key relations in equilibrium 

Let us consider site i and sites j G B { in a random graph Q G G(c, N). Here, we assume 
that the sites j, respectively, can be regarded as the root of the Cayley tree T*, of which 
generation is not increased to the direction of site i. 

Under this condition, when we ignore site % for site j G Bi, we can regard 
Zj '(er)/ ^2 a , Zj(a') as the probability that <7j takes a value a. Here, let us restate the 
bonded sites of site i, {ji, ■ ■ ■ ,ji t }, and the unbonded sites of site i, {j^+i, ■ ■ ■ ,j c }- Using 
this representation, we can obtain the concrete expression of the probability Pf ,cq (a,u) 
that with a realization of J, (<7j,Wj) takes (er, u) in equilibrium as follows: 

11 v Z J (o-') J-l v Z J ((t'Y K ' 

where Nq is determined by J2au Pf' eq (a, u) = 1. 
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Let us consider the probability p° q u defined as 



p3,« = *(£ J «-0^ ,e >>«)> ( R2 ) 

jeBi 

which corresponds to p a ,i,u(t > *eq)- Here, we have assumed that p° q u does not depend 
on the chosen site i, and Zj(—1)/Zj(1) does not depend on site j. When we define 
exp(/32/i av ) = Zj(—1)/Zj(l), using Q(h) which has been already obtained in (A. 3), 
we can find exp(/32/i av ) = J dhQ(h) exp(/32h). Then, we obtain the following exact 
expression of the probability p c ^ lu . 

^^(0 p,(i - pr '(i, 

exp(-/3cr(7 - 2m)) exp(w/32/i av )(l + exp(/32/i av )) c ~\ (B.3) 

where we used relation (12) and determined N v by the condition E u <z P^i u = 1- 

Next, we discuss the joint probability Pfj' eq ((o~, u), (a', u')) that in equilibrium, with 
a realization of J, (cr^Uj) and (aj,Uj) take (a,u) and (a',u'), respectively Let us 
consider the probability 

iff(foi, M ),(y,iV))= 

tf(J - - 1)*( J2 J ™ - l ) 5 ( E J *» " ^' eq (M, (<",*')), (B.4) 

which corresponds to Pb((o~,1,u), (a' ,1' ,u');t > t eq ) assuming P^((a,l,u),(a',l',u')) 
does not depend chosen bonded sites In the same way as p° q u , we obtain the exact 
expression as follows: 

^ - ( „< -V + 1) ) ( „ - 1) ) "^(-^^ 

exp(-/V(Z' - 1 - 2(u' - 5(a + 1))) - pa (I - 1 - 2(u - + 1))) 

+(u' - (J(o- + 1) + u - + l))/32/i av )(l + exp(/32/i av )) 2c ^- r , (B.5) 

where N e is determined by J2a E„<z E U '</' ^bKK, u '), M )) = L It should 
be noted that from expressions (B.3) and (B.5), we derive key relations 

P°\(+l,l',u'),(a,l,u)) _ (l-u)p c ? u 



E,» „» ^ q ((+i, J', «0, «")) E,» „«(i" - «>3» „* ' 

iff((-l,Z',n')>,i,rQ) _ upl\ u 

E, v « ^ q ((-i, «'), (*> E,» «» «^3„ „„ ' 
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